function fig7_21

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

figure(1)
clf

% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
hold on

tmax=500;
y0=[0.1 0 2 2];

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);
%[t,y] = ode23s(@rhs,[0 tmax],y0);
z=y(:,4)-y(:,2);

[tt,yy] = ode45(@rhs2,[0 tmax],y0);
%[t,y] = ode23s(@rhs,[0 tmax],y0);
zz=yy(:,4)-yy(:,2);

plot(t,z,'b','LineWidth',1.2)
plot(tt,zz,'r','LineWidth',1.2)

axis([0 500 -0.2 4])
% commands to label each axes

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$t$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('$\phi$\,-\,axis','Interpreter','latex','FontSize',14)

set(gca,'xtick',[0 100 200 300 400 500])

grid on
legend({'$\, a=1$','$\, a=-1$'},'Location','East','Interpreter','Latex','FontSize',12,'FontWeight','bold','AutoUpdate','off')

box on
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/fig620x.pdf')

return

figure(2)
clf
set(gcf,'Position', [993 1123 666 247])

hold on
plot(t,y(:,1),'k')
plot(t,y(:,3),'--k')

%plot(tt,yy(:,1),'k')
%plot(tt,yy(:,3),'--k')

axis([0 10 0 2]);
% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('r - axis','FontSize',14,'FontWeight','bold')

%set(gca,'xtick',[0 100 200 300 400 500]);

grid on
% command to put legend into plot
legend(' r_1',' r_2','Location','NorthEast');

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14);




%  define f(t,y)
function dy=rhs(t,y)
global a0
% y(1)=r(1) y(2)=theta(1)  y(3)=r(2)  y(4)=theta(2)
dy=zeros(4,1);
ep=0.01;
a=1;
b=1;
a1=1; b1=1;  k1=1;
a2=a1;  b2=b1;  k2=k1;

z1=y(3)*cos(y(4)) - y(1)*cos(y(2));
z2=y(3)*sin(y(4)) - y(1)*sin(y(2));

f1 = -a*z1 - b*z2;  g1 = b*z1-a*z2;
f2 = a*z1 + b*z2; g2=-b*z1+a*z2;

F1=f1*cos(y(2))+g1*sin(y(2));
G1=(g1*cos(y(2))-f1*sin(y(2)))/y(1);

F2=f2*cos(y(4))+g2*sin(y(4));
G2=(g2*cos(y(4))-f2*sin(y(4)))/y(3);

dy(1)=-b1*y(1)*(y(1)^2-k1^2)+ep*F1;
dy(2)=-a1+ep*G1;
dy(3)=-b2*y(3)*(y(3)^2-k2^2)+ep*F2;
dy(4)=-a2+ep*G2;



%  define f(t,y)
function dy=rhs2(t,y)
global a0
% y(1)=r(1) y(2)=theta(1)  y(3)=r(2)  y(4)=theta(2)
dy=zeros(4,1);
ep=0.01;
a=-1;
b=1;
a1=1; b1=1;  k1=1;
a2=a1;  b2=b1;  k2=k1;

z1=y(3)*cos(y(4)) - y(1)*cos(y(2));
z2=y(3)*sin(y(4)) - y(1)*sin(y(2));

f1 = -a*z1 - b*z2;  g1 = b*z1-a*z2;
f2 = a*z1 + b*z2; g2=-b*z1+a*z2;

F1=f1*cos(y(2))+g1*sin(y(2));
G1=(g1*cos(y(2))-f1*sin(y(2)))/y(1);

F2=f2*cos(y(4))+g2*sin(y(4));
G2=(g2*cos(y(4))-f2*sin(y(4)))/y(3);

dy(1)=-b1*y(1)*(y(1)^2-k1^2)+ep*F1;
dy(2)=-a1+ep*G1;
dy(3)=-b2*y(3)*(y(3)^2-k2^2)+ep*F2;
dy(4)=-a2+ep*G2;

